library(rio)
library(ggplot2)

data <- import("estimation_gamma_vars.dta")

if (!("gamma" %in% names(data))) data$gamma <- data$experience

data <- subset(data, select = c(starty, endy, ruler_unique, gamma, polity2))

data <- do.call("rbind", (apply(data, 1, function(x) {
  
  data.frame(
    year = seq(x["starty"], x["endy"]),
    ruler = x["ruler_unique"],
    gamma = x["gamma"],
    polity = x["polity2"]
  )
})))

data <- subset(data, year < 2018 & year >= 1950)

data$democracy <- ifelse(as.numeric(as.character(data$polity)) >= 6, "Democracy", "Non-democracy")

data$gamma <- as.numeric(as.character(data$gamma))

annualMean <- aggregate(gamma ~ year, data, mean)
annualQuantileLow <- aggregate(gamma ~ year, data, quantile, .2)
names(annualQuantileLow)[2] <- "low"
annualQuantileHigh <- aggregate(gamma ~ year, data, quantile, .8)
names(annualQuantileHigh)[2] <- "high"

annual <- merge(annualMean, annualQuantileLow, by = "year")
annual <- merge(annual, annualQuantileHigh, by = "year")

print(p1 <- ggplot(annual, aes(x = gamma, y = year)) +
  geom_errorbarh(aes(xmin = low, xmax = high), height = .5, col = "gray") +
  geom_point(size = .8) +
  theme_minimal() +
  labs(x = "Average PolEx", y = "Year") +
  coord_flip())

annualMean <- aggregate(gamma ~ year + democracy, data, mean)
annualQuantileLow <- aggregate(gamma ~ year + democracy, data, quantile, .2)
names(annualQuantileLow)[3] <- "low"
annualQuantileHigh <- aggregate(gamma ~ year + democracy, data, quantile, .8)
names(annualQuantileHigh)[3] <- "high"

annual <- merge(annualMean, annualQuantileLow, by = c("year", "democracy"))
annual <- merge(annual, annualQuantileHigh, by = c("year", "democracy"))

dodge <- position_dodge(width = .8)
print(p2 <- ggplot(annual, aes(y = gamma, x = year, col = democracy)) +
        geom_errorbar(aes(ymin = low, ymax = high), width = .5, position = dodge) +
        geom_line(size = 1.2) + # geom_point(size = 1.2, position = dodge) +
        theme_minimal() +
        labs(y = "Average PolEx", x = "Year") +
        guides(col = guide_legend(title = "")) +
        scale_color_manual(values = c("darkgray", "black")) +
        theme(legend.position = "bottom") )

print(p3 <- ggplot(annual, aes(y = gamma, x = year, col = democracy)) +
        geom_line(size = 1.2) +
        theme_minimal() +
        labs(y = "Average PolEx", x = "Year") +
        guides(col = guide_legend(title = "")) +
        scale_color_manual(values = c("darkgray", "black")) +
        theme(legend.position = "bottom") )

cairo_pdf("../manuscript/figures/average_gamma_over_time.pdf", width = 7, height = 4, onefile = TRUE)
print(p1)
print(p2)
print(p3)
dev.off()

